function [results, standard_deviation, img_points] = read_circ_target(filepath, varargin)
% READ_CIRC_TARGET, a MATLAB target measuring script.
% Copyright (C) 2012 Glenn D. Sweeney
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.
% 
%
%
%
% Description:
% READ_CIRC_TARGET is a function designed for measuring average digital
% counts for each patch in an image of a circular target. It uses user
% input to locate the target in the image, and is capable of processing
% oblique images of targets.
%
% Usage:
% READ_CIRC_TARGET(filepath) opens a file or folder specified by 'filepath'
% and produces a GUI for manual target identification. The user must click
% in the center areas between patches at each of the four cardinal
% directions in the following order: North, East, South, West. (see
% illustration below for an example of correct position of click)
%
%                          ___________________
%                          |        |        |
%                          | (top   |  (top  |
%                          |  left) *  right)|
%                          |        |        |
%                          |________|________|
%
% 
% Mandatory Argument:
% 
%   filepath            This argument specifies the location of the image
%                       (or folder of images, see below) to be processed
%
% Optional Arguments:
%   Optional arguments are supplied using a 'parameter', (value) syntax,
%   where function parameters may be specified in any order, and then the
%   paired value (if any, depending on the parameter) is passed in as the
%   next function argument. For example, if the parameter 'CropFactor' is
%   to be set to 0.5, the calling syntax would be:
%
%   READ_CIRC_TARGET(filepath, 'CropFactor', 0.5)
%
%   A table of parametrs with a description of acceptable values is
%   provided below:
% _________________________________________________________________________
% |                    |                                                  |
% |    'CropFactor'    | The value immediately following 'CropFactor'     |
% |                    | should be a decimal between zero and one. This   |
% |                    | value denotes the proportion  of each patch to   |
% |                    | be removed from calculation.                     |
% |                    | Default value is 0.40.                           |
% |                    |                                                  |
% |     'Mosaic'       | Mosaic allows a user to input a 2x2 integer      |
% |                    | array denoting the mosaic pattern present in a   |
% |                    | grayscale image. The pattern must consist only   |
% |                    | of 1, 2, and 3, corresponding to R,G,B channels. |
% |                    | An example is as follows:                        |
% |                    |      bayerarr =    ['2' '1'; ...                 |
% |                    |                     '3' '2']                     |
% |                    |      read_circ_target(...,'Mosaic', bayerarr)    |
% |                    | If 'Mosaic' is not specified, it is assumed that |
% |                    | the image has already been demosaiced.           |
% |                    |                                                  |
% |    'Normalize'     | Normalize denotes the method of graphical        |
% |                    | display used. It does not affect output data.    |
% |                    | This tag is NOT followed by a value.             |
% |                    |                                                  |
% |   'NudgeFactor'    | Setting this will artifically expand or contract |
% |                    | the measurement circle to match to centers.      |
% |                    | Default value is 1.05.                           |
% |                    |                                                  |
% |    'NumPatches'    | Change the number of patches to be measured.     |
% |                    | Default is 12.                                   |
% |                    |                                                  |
% |    'PatchRatio'    | Defines the size of each patch as a fraction of  |
% |                    | the diameter of the circle being measured.       |
% |                    | Default (standard for 12-patch) is 0.2.          |
% |                    |                                                  |
% |     'Folder'       | Specifying the folder parameter will process a   |
% |                    | folder of identically-composed images. This tag  |
% |                    | is not followed by a value.                      |
% |                    | Usage: read_rect_target(folderpath, [r,c], ...,  |
% |                    |          'Folder')                               |
% |                    |                                                  |
% |    'ImgPoints'     | Allows the user to pass in a set of predetermined|
% |                    | control points for processing (same format as    |
% |                    | imgpoints return value)                          |
% |                    |                                                  |
% |  'SuppressDraw'    | When this flag is set, no figures will be created|
% |                    | for output .This parameter is not followed by a  |
% |                    | value.                                           |
% |____________________|__________________________________________________|
%
% Output Values:
%   READ_CIRC_TARGET returns three sets of values. The first is an RxC
%   array, where each array location has a vector of average digital counts
%   for each channel. Accordingly, for an RGB image, this output will be a
%   3d array, indexable by: (row, col, chan).
%   Results will always be doubles, but scaled to match the bit depth of
%   the input image.
%   Usage: result = READ_CIRC_TARGET(...)
%
%   Optionally, READ_CIRC_TARGET will return standard deviations for each
%   patch, in the same format as patch digital counts.
%   Usage: [digital_count, stdev, ~] = READ_CIRC_TARGET(...)
%
%   Finally, if many images are going to be processed with the same target,
%   but for some reason the 'Folder' parameter is unsatisfactory, the
%   selected coordinates are also returned, and can be passed into a
%   subsequent call of READ_RECT_TARGET through the 'ImgPoints' parameter.
%   Usage: [digital_count, stdev, impoints] = READ_CIRC_TARGET(...)




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%  Set Up Function Parameters  %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Default variable states
crop_factor = .25;
mosaic_array = 0;
peak_normalize = 0;
num_patches = 12;
folder_flag = 0;
points_flag = 0;
draw_flag = 1;
ratio_set_flag = 0;
nudge_set_flag = 0;


% Check all arguments
i = 1;
while i <=numel(varargin)
    
    switch varargin{i}
        case 'CropFactor'
            crop_factor = varargin{i+1};
        case 'Mosaic'
            mosaic_array = varargin{i+1};
        case 'Normalize'
            peak_normalize = 1;
            % Don't need the second input, so back off by 1
            i = i-1;
        case 'NudgeFactor'
            nudge_factor = varargin{i+1};
            nudge_set_flag = 1;
        case 'NumPatches'
            num_patches = varargin{i+1};
        case 'PatchRatio'
            patch_ratio = varargin{i+1};
            ratio_set_flag = 1;
        case 'Folder'
            folder_flag = 1;
            i = i-1;
        case 'ImgPoints'
            points_flag = 1;
            img_points = varargin{i+1};
        case 'SuppressDraw'
            draw_flag = 0;
            i = i-1;
        otherwise
            error('Unknown function argument.');
    end
    
    i = i+2;
end

% Preset ratio depending on number of patches
if ratio_set_flag == 0
    if num_patches == 20
        patch_ratio = 0.10;
    else
        patch_ratio = 0.2;
    end
end

% Preset nudge factor depending on number of patches
if nudge_set_flag == 0
    if num_patches == 20
        nudge_factor = 1.02;
    else
        nudge_factor = 1.05;
    end
end

% Get all files at the specified directory location
file_list = dir(filepath);

% on a mac, hidden files begin with '.', and should be ignored
while strcmp(file_list(1).name(1), '.')
    file_list = file_list(2:end);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%  Read and Process First Image  %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% Read in the first target image
if folder_flag
    input = imread([filepath file_list(1).name]);
else
    input = imread(filepath);
end


% Process the mosaic by splitting across color channels. Populate all other
% pixels with -1 to avoid computational biasing.
if mosaic_array == 0
    target = input;
else
    tar_size = size(input);
    tar_height = tar_size(1);
    tar_width = tar_size(2);
    target = zeros([tar_size, 3])-1;
    target(1:2:tar_height, 1:2:tar_width, mosaic_array(1,1)) = input(1:2:tar_height, 1:2:tar_width);
    target(2:2:tar_height, 2:2:tar_width, mosaic_array(2,2)) = input(2:2:tar_height, 2:2:tar_width);
    target(1:2:tar_height, 2:2:tar_width, mosaic_array(1,2)) = input(1:2:tar_height, 2:2:tar_width);
    target(2:2:tar_height, 1:2:tar_width, mosaic_array(2,1)) = input(2:2:tar_height, 1:2:tar_width);
end


% Calculate the size of the target
tar_size = size(target);
tar_height = tar_size(1);
tar_width = tar_size(2);

% Test for number of color channels
if numel(tar_size) == 3
    num_channels = tar_size(3);
else
    num_channels = 1;
end


% Define the function output
output = zeros(num_patches, num_channels);
stdev = zeros(num_patches, num_channels);

results = zeros(num_patches, num_channels, numel(file_list));
standard_deviation = zeros(num_patches, num_channels, numel(file_list));


if ~points_flag
    
% Display the image, and read in 4 corners specified by user
if peak_normalize == 0
    imshow(input);
else
    imshow(input, []);
end
hold on
num_good_points = 0;
    
    % first row: x   Second row: y    3rd row: homogeneous coordinate
    img_points = zeros(3, 4);
    
    % Size of the corner dots to be used (ceil to ensure at least one pixel)
    % Defined as 1/100 the smaller dimension
    dot_size = ceil(min([tar_height tar_width])/100);
    
    % Get target registration points
    while num_good_points < 4
        
        [x_click y_click button] = ginput(1);
        % If it was left mouse click:
        if button == 1
            % Show the pretty red squares
            fill([x_click+dot_size x_click+dot_size x_click-dot_size x_click-dot_size], ...
                [y_click+dot_size y_click-dot_size y_click-dot_size y_click+dot_size], 'r');
            num_good_points = num_good_points + 1;
            img_points(1, num_good_points) = x_click;
            img_points(2, num_good_points) = y_click;
            img_points(3, num_good_points) = 1; % Generating homogeneous coordinates
            drawnow;
        end
    end
    
    % Close the current window
    close all;
    
end

% Sorts the four image points so they go clockwise from top
% FOR NOW, ASSUME THAT IS HOW THEY ARE SELECTED.

% Target points assume upper left corner is zero, and bottom right is
% (1,1), clockwise from top center
target_points = [0  1  0 -1; ...
    1  0 -1  0; ...
    1  1  1  1];

% Generate A matrix for each correspondance:
% Projecting from the src (source) system into the dst (destination) system
A = [];
for i = 1:4
    src = target_points(:,i);
    dst = img_points(:,i);
    
    A = vertcat( A, [0,0,0,   -1*src',   dst(2)*(src');...
        src',    0,0,0,   -1*dst(1)*(src')]);
    
end

% Obtain SVD of A
[~,~,V] = svd(A);

% Find H using last column of V (the null space of A)
transform_mat = reshape(V(:,end), 3, 3)';
transform_mat = transform_mat./transform_mat(9);

patch_angle = 360/num_patches;

% Generate sample coordinates for each angular position
sample_coords = [sin(((patch_angle/2):patch_angle:360)/180*pi);cos(((patch_angle/2):patch_angle:360)/180*pi)].*nudge_factor;






% Iterate through each image
for imnum = 1:numel(file_list)
    
    
if folder_flag
    input = imread([filepath file_list(imnum).name]);
else
    input = imread(filepath);
end
    
    
    
    % Guess at the bit depth of the image.
    if isa(input, 'uint8')
        tar_bitval = 8;
        input = im2double(input);
    elseif isa(input, 'uint16')
        tar_bitval = 16;
        input = im2double(input);
    else
        tar_bitval = ceil(log(double(max(input(:))))/log(2));
        input = input./(2^tar_bitval-1);
    end
    
    
    % Process the mosaic by splitting across color channels. Populate all other
    % pixels with -1 to avoid computational biasing.
    if mosaic_array == 0
        target = input;
    else
        tar_size = size(input);
        tar_height = tar_size(1);
        tar_width = tar_size(2);
        target = zeros([tar_size, 3])-1;
        target(1:2:tar_height, 1:2:tar_width, mosaic_array(1,1)) = input(1:2:tar_height, 1:2:tar_width);
        target(2:2:tar_height, 2:2:tar_width, mosaic_array(2,2)) = input(2:2:tar_height, 2:2:tar_width);
        target(1:2:tar_height, 2:2:tar_width, mosaic_array(1,2)) = input(1:2:tar_height, 2:2:tar_width);
        target(2:2:tar_height, 1:2:tar_width, mosaic_array(2,1)) = input(2:2:tar_height, 1:2:tar_width);
    end
    
    if draw_flag
    
        if peak_normalize == 0
            imshow(target);
        else
            imshow(target./max(target(:)));
        end

        hold on
    end
    
    % Iterate through each patch of the target
    for sample = 1:num_patches
        x_patch = sample_coords(1,sample);
        y_patch = sample_coords(2, sample);
        
        % Define the current patch in the normalized coordinate system
        target_coords = [(x_patch - patch_ratio*(1-crop_factor)), (x_patch + patch_ratio*(1-crop_factor)), (x_patch + patch_ratio*(1-crop_factor)), (x_patch - patch_ratio*(1-crop_factor)); ...
            (y_patch + patch_ratio*(1-crop_factor)), (y_patch + patch_ratio*(1-crop_factor)), (y_patch - patch_ratio*(1-crop_factor)), (y_patch - patch_ratio*(1-crop_factor)); ...
            1                    ,                       1                ,                      1                 ,                      1                ];
        
        % Transform to target coordinate system
        img_coords = transform_mat*target_coords;
        
        % Convert from homogeneous coordinates to the images
        patch_coords = [img_coords(1,:)./img_coords(3,:); img_coords(2,:)./img_coords(3,:)];
        
        if draw_flag
        
            % Draw the pretty red boxes
            line = plot([patch_coords(1,:), patch_coords(1,1)], [patch_coords(2,:), patch_coords(2,1)], '-r', 'LineWidth', 2);

            % Force a re-draw
            drawnow;
        end
        
        % Create a ROI mask
        mask = roipoly(target, patch_coords(1,:), patch_coords(2,:));
        
        for i = 1:num_channels
            
            % For each color channel, apply the mask
            patch = target(:,:,i);
            patch = patch(mask);
            patch = patch(patch ~= -1);
            
            % Mean the ROI and report
            output(sample, i) = mean(patch(:));
            stdev(sample, i) = std(patch(:));
            
        end
        
    end
    
    % Re-order patches properly
    % Fold in half in the center, then reshape
    
    outputOrder = [rot90(1:(num_patches/2)), rot90((num_patches/2+1):(num_patches), -1)];
    outputOrder = outputOrder';
    outputOrder = outputOrder(:);
    
    output = output(outputOrder, :);
    stdev = stdev(outputOrder, :);
    
    % Scale output patches to original bit depth
    output = output.*(2^tar_bitval-1);
    stdev = stdev.*(2^tar_bitval-1);
    
    results(:,:,imnum) = output;
    standard_deviation(:,:,imnum) = stdev;
    
end

end

